September 23 2021
The lab counts for up to 4 points towards the final grade of the course.
Read the instructions below step by step. Copy and paste each chunk of code at a time in your R script. Select the code with your mouse or shift/arrow keys and then run it by pressing the keys control-enter simultaneously. Complete the lab deliverables document and upload it to canvas along with the requested pdf file.
Have fun!
Reproduce the steps for data downloading described in lab 2 to download from Earth explorer (http://earthexplorer.usgs.gov), the Landsat image that covers the extent of the island of Guam, US acquired on 2021/02/08 (path 100, Row 051).
Open R studio, and load required libraries.
Setup the folder where you stored the downloaded image as your working directory. Make sure you use forward slashes in the wd name.
library(raster)
library(rgdal)
library(RStoolbox)
wd="/Users/tug61163/Documents/Courses/IntroRemoteSensing/2021Fall/Class4/Data"
setwd(wd)
dir(wd)
Check the names of the files in your working directory. Then copy the name of the compressed (.tar) Landsat file name and paste it in the untar function to decompress the file.
The list.files() function produces a file with the names of all the files in the working directory that share a common pattern. In this case, it will list all files that end with the extension .TIF.
Identify the position of the bands of interest in the tifs object and enter them in the tifBandNames.Then use them to stack the bands and assign a name to the bands within the stack:
untar("LC08_L2SR_100051_20200208_20201016_02_T1.tar")
tifs <- list.files('.', pattern='.TIF')
tifs # lists all the files in the working folder with the extension .TIF
tifBandNames=tifs[c(3:9,1)] # enter the index position for the bands of interest
L8=stack(tifBandNames)
names(L8)=tifBandNames
names(L8)
plotRGB(L8, r=5, g=4, b=3, stretch="lin")
e=drawExtent()
L8rsz=crop(L8, e)
plotRGB(L8rsz, r=5, g=4, b=3, stretch="lin")
## To read only QA classes, use:
plot(L8rsz[[8]])
click(L8rsz[[8]], n=5, cell=TRUE)
The numbers assigned to different pixels in the QA layer represent different quality attributes defined by the bit values that reproduce such number. You can see the interpetation in Tables 6-2 and 6-3 (pages 13 and 14) here: https://www.usgs.gov/media/files/landsat-8-collection-2-level-2-science-product-guide.
As an example, let’s decode the value of a pixel representing 1. cloud, 2. cloud shadow, 3. clear land.
Before that, I have created a function to transform a number in a binary representation in the format used by the Landsat Science Team. Let’s run it to load it to the environment
int2bits=function (x, n, nbits=16)
# Converts an integer into the corresponding bit denomination for a number of bits specified by n
# x: integer to be converted into bits
# n: bit numbers to be retreived. It can be a number or a vector specifying the bit positions to retrieve.
# The first bit is counted as zero and they are read from right to left
{
bit <-intToBits(x)
rev(as.numeric(paste(tail(rev(as.integer(bit)), nbits))))[n+1]
#if (reverse==TRUE){bit=rev(bit)}
#return(bit)
}
Let’s decompose and interpret the QA information embedded in some selected numbers (you can explore other numbers that you obtained when you used the click() function:
int2bits(22280, c(0:15)) # Cloud
int2bits(23888, c(0:15)) # Cloud shadow
int2bits(21824, c(0:15)) # Clear observation of land
int2bits(21888, c(0:15)) # Clear observation of water
We will use those bit values to create masks representing clouds and cloud shadows. The package RStoolbox has a function to classify different components of the QA layer (classifyQA()), unfortunately it only works for Landsat collection 1 which is going to be unavailable very soon. Let’s instead, load a function that I created for that purpose:
bitMask=function(r, bit){
# produces a mask that excludes as NA, pixels with a value of 1 in the bit position designated in the argument "bit" in the landsat QA layer. All other pixels are assigned a value of 1
x <- calc(r, fun = function(x) int2bits(x, bit))
v <- raster::getValues(x)
v[v==1]=NA
v[v==0]=1
x <- raster::setValues(x, v)
return(x)
}
The bitMask function creates a mask that excludes pixels in the QA layer with a value of 1 in a bit position designated by the argument “bit”. For instance, pixel values with a value equal to 1 in bit 3 (correspond to clouds) will be re-classified as NA and all other pixels will be reclassified as 1. This function might take several minutes to process so be patient.
clouds=bitMask(L8rsz[[8]],3)
plot(clouds)
Multiplying the raster image by the cloud mask, removes areas covered by clouds.
L8rszMskd=L8rsz*clouds
plotRGB(L8rszMskd, r=5, g=4, b=3, stretch="lin")
Let’s remove cloud shadows. For that purpose we use the fourth bit instead of the third.
shadows=bitMask(L8rsz[[8]],4)
plot(shadows)
Multiplying the raster image by the shadow mask, removes shadowed areas.
L8rszMskd2=L8rszMskd*shadows
plotRGB(L8rszMskd2, r=5, g=4, b=3, stretch="lin")
pdf("VGutierrez_Lab3.pdf")
plotRGB(L8rsz, r=5, g=4, b=3, stretch="lin")
plotRGB(L8rszMskd, r=5, g=4, b=3, stretch="lin")
plotRGB(L8rszMskd2, r=5, g=4, b=3, stretch="lin")
dev.off()
Define your study area, including the full island similar to this extent.
Produce a single pdf file with three RGB images as shown in step 5 of the lab: 1. Unmasked image, 2. image with clouds masked out, 3. Image with clouds and cloud shadows masked out. Upload the pdf to module 4 in canvas (3 pts).
For the numbers below, mark with an x the corresponding classes (some values can have more than one) in the QA layer produced for Landsat 8, collection 2, level 2. Upload your answers to module 4 in canvas (1 pt):
21952 ___ dilated cloud, ___ cirrus, ___ Cloud, ___ Shadow,___ Water
22018 ___ dilated cloud, ___ cirrus, ___ Cloud, ___ Shadow,___ Water
23826 ___ dilated cloud, ___ cirrus, ___ Cloud, ___ Shadow,___ Water
55052 ___ dilated cloud, ___ cirrus, ___ Cloud, ___ Shadow,___ Water